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Turbulence modeling for non-equilibrium flows 


By P. A, Durbin 


1. Motivation 

Two projects are reported herein. The first is the development and testing of 
an eddy viscosity transport model. This project also is a starting point for our 
work on developing computational tools for solving turbulence models in complex 
geometries— the computational aspect is collaborative with Nagi Mansour. The 
second project is a stochastic analysis of the realizability of Reynolds stress trans- 
port models. This work was motivated by private discussions with Charles Speziale 
(Boston University) and related collaborative work with him, and also by Steve 
Pope’s presentation at the (Bill) Reynolds Turbulence Symposium in Monterey this 
year (Pope 1993). 


1.1 Eddy viscosity transport 

Momentum mixing-length, eddy viscosity models have been popular since the idea 
was introduced by Prandtl. They invoke a quasi-equilibrium assumption, essentially 
requiring that the mixing of lumps of fluid takes place in a time short compared 
to that of the mean flow evolution. In many situations of fluid dynamical interest, 
mixing cannot be assumed instantaneous and a dynamical equation for the turbu- 
lence is required. A rather important example is that of boundary layers subjected 
to strong pressure gradients: on a practical level, it has been found that mixing 
length models are unable to predict this type of strongly non-equilibrium flow. A 
reasonable step is to formulate an analytical eddy viscosity transport model that 
extends the mixing- length idea by admitting non-equilibrium effects; this is what I 
have done. 

This research has a strong practical incentive. Engineering fluid dynamicists 
are recognizing a need for more elaborate turbulence models as the complexity of 
the flows they calculate increases. This has led to a willingness, even a desire, 
to introduce turbulent transport models into prediction codes. Currently most 
aerodynamics codes (that solve the Reynolds averaged Navier-Stokes equations) 
use algebraic eddy viscosity models. It has been found that algebraic models are 
unable to predict the region of reversed flow near the trailing edge of an airfoil at 
angle of attack. The failure in this particular case indicates a general inadequacy for 
calculating complicated mean flows. This shortcoming of simple algebraic models 
motivated recent research into dynamical equations for eddy viscosity, including that 
reported herein (Durbin et ai 1994). Baldwin and Barth (1990) first proposed that 
an eddy-viscosity transport equation might be effective in complex aerodynamic 
flows. The Baldwin-Barth study led to further development of their formulation by 
Spalart and Almaras (1992) and was a primary impetus for the present work. 
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Notable features of the present model are that it uses an elliptic relaxation equa- 
tion to avoid damping functions, that it is formulated solely in terms of local vari- 
ables, and that it is tensorially and Galilean invariant. In these respects, the model 
was formulated with complex flows in mind. 

1.2 Realizability 

The exact, unclosed Reynolds- stress transport equations are usually the start- 
ing point for formulations of second-moment closure models. Modeling consists 
of replacing unclosed terms by semi-empirical formulae that express these terms 
as functions _of the dependent variables. After introducing models, quantities with 
names like 4 u 2? no longer represent non-negative functions obtained by squaring and 
averaging a random variable; rather, they are simply the dependent variables of the 
model; they are obtained as the solution to a differential equation. However, it is 
desirable to formulate the equations of the model so that variables like u 2 do main- 
tain their non- negativity. In essence, this is the issue of realizability in second-order 
turbulence closure modeling. 

The present report demonstrates how realizability can be addressed by a con- 
structive method. This involves formulating a stochastic process for which the 
Reynolds-stress model is the exact evolution equation of second moments. The 
model then is guaranteed to be realizable because it is exact for a well defined 
stochastic process. In the present analysis, second-moment closure models of the 
type currently in use are shown to be exact for the statistics of a particular form 
of Langevin equation. When that Langevin equation is well-defined, the Reynolds- 
stress model is guaranteed to be realizable. 

It is assumed that the appropriate physics are accommodated by the moment 
closure; the stochastic analysis is purely a mathematical method for analyzing such 
models. The realizability criteria that are derived are sufficient , but not necessary, 
conditions; also, the analysis is only of homogeneous turbulence. 

2. Accomplishments 

2.1 The eddy viscosity transport model 

The model is described at length in Durbin et al. (1994). It consists of a parabolic 
transport equation for the eddy viscosity: 

D v 2 

— v T = V(l/ + V'f)VvT + Pv “ C 2 \S\ut - C 4 ^2 ; (1) 

and an elliptic relaxation equation for P v : 

L 2 p V 2 P v -P v = c z \Vvt \ 2 -\S\v T . (2) 

The dependent variable P u contains the turbulence production. The elliptic relax- 
ation in (2) introduces a wall effect that suppresses the production of eddy viscosity 
near a surface. In (1) and (2) the c's are model constants, the Vs are length scales, 
and |S| 2 = 1 1 2 ( d % Uj + djUi)(djUi + d{Uj) is a measure of the mean rate of strain. 
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Figure 1. Contours of constant U, showing trailing edge separation and flow 
geometry. 

The ellipticity in (2) is important near surfaces; away from surfaces, or regions of 
strong inhomogeneity, the model relaxes to a parabolic transport equation. The first 
term on the right side of (1) describes turbulent and molecular transport; the third 
and fourth terms model dissipation of u T , as does the first term on the right side 
of (2); the second term on the right side of (2) allows for production of turbulence 
from mean flow gradients. 

To complete the model, the length scales must be prescribed and then the con- 
stant coefficients chosen. A difficulty of a model with only one primary dependent 
variable— or, more precisely, only one non-negative dependent variable— is that lo- 
cal turbulent length and time scales cannot be formed from turbulence quantities 
alone. For this reason L p and L„ have been prescribed as functions of the local 
mean rate of strain, as well as of vt'- 

r 2 |S| 2 . IVt/rl 2 

L * - |VS| 2 +Cm_ [sp~ (3) 

L p = c pmin(l£,max(i/r,c 2 i/)/|S|) . ( 4 ) 

The boundary conditions to the model axe 

vt = n ■ Vut = 0 (5) 

at a no-slip surface with unit normal n, 

h ■ Vu T = P v = 0 (6) 

at the edge of a boundary-layer, or prescribed free-stream values far from an airfoil. 
No surface boundary condition needs to be imposed on P„ to satisfy (5). 
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FIGURE 2. Mean velocity profiles at locations along the suction surface of an 
airfoil at angle of attack, showing trailing edge separation. The profiles are at 
x/c = 0.62, 0.675, 0.731, 0.786, 0.842, 0.897 and 0.953. The circles are experimental 
data. 

The selection of model constants is described in Durbin et al. (1994). To sum- 
marize, the present values are C 2 — 0.85, c 4 — 0.2, cj = 3.3, c p = 1.2 c m = 2 and 
c 3 = (1 - c 2 )//c 2 + 1 - c 4 /( 1 + K*c m ) where k is the von Karman constant. 

Figures 1 and 2 illustrate the model solutions. They are from a computation 
of flow around a NACA4412 airfoil at 14° angle of attack. Contours of constant 
t/-component of mean velocity are plotted in figure 1. These illustrate that the flow 
separates just upstream of the trailing edge. The whole computational domain is 
not shown: the inflow is 14 chords upstream and is angled upward by 14°. 

Figure 2 contains mean velocity profiles on the suction surface near the trail- 
ing edge. The last two profiles show that the model predicts the correct degree of 
separation. The Baldwin-Lomax, equilibrium eddy-viscosity model fails to predict 
separation in this flow; the rapid evolution of the flow as it undergoes separation 
violates the assumption of near equilibrium. The operation of the non-equilibrium 
model would seem to be as follows: the eddy viscosity transports mean momentum 
from the free-stream to the surface, counter balancing deceleration by the adverse 
pressure gradient. As the boundary layer thickens, eddy viscosity will be produced 
in the outer region of the boundary layer. However, this process is not instan- 
taneous. If the pressure gradient decelerates the near-wall flow before the eddy 
viscosity can increase, the flow will separate. The transport equation describes 
the temporal evolution and spatial redistribution of the eddy viscosity as the flow 
passes through separation. Equilibrium, or algebraic, models assume that the eddy 
viscosity adjusts instantaneously to the level of shear and hence they overpredict 
the turbulent transport in this trailing edge flow. 
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2.2 Stochastic analysis of realizability 

The stochastic method will be presented by describing how it is used to derive a 
sufficient condition for realizability of the most general, linear, second-order closure 
model for Reynolds stresses. However, this method is far more widely useful and can 
be applied to all existing Reynolds stress transport models, including those that are 
non-linear in the Reynolds stress tensor. The method is presented at greater length 
in Durbin and Speziale (1993). It involves showing the Reynolds-stress model to be 
the exact second-moment equation of a Langevin stochastic differential equation. 
The exact moment equation of a well defined stochastic process is realizable, by 
definition — for instance, quantities like u 2 have non-negative values that could be 
obtained by squaring and averaging a random variable. 

In homogeneous turbulence, the General Linear Model consists of the following 
ordinary differential equation for the evolution of the Reynolds stress tensor 

- |m#) ~ C *( P H ~ \ P6 H) 

~ c 3 {Dij — — C»kSij + Pij — -fijj. (7) 

In this equation, P,j = — UiU k d k Uj — UjU k d k Ui is the production tensor, D t j = 
—UiUkdjUk — UjUkdiUk is a tensor introduced by Launder et al. (1975), P = x j 2 P x i 
is the rate of turbulent energy production, and e is its rate of dissipation. All the 
Cj’s are empirical, numerical constants and T is a time-scale. 

We will show (7) to be the second moment of the stochastic differential equation 


dUi = ~2 T u ' di + ( C2 “ 1 hkdkUidt + c 3 u k diU k dt 

+ y/c^ed\Vi(t) + y/c7eMi k dW' k (i) (8) 


with 




— 1 + (C2 + C3 )“- 



(9) 


In (8), W and W' are independent Weiner processes (dW^dVV; = 0). M is a 
symmetric matrix that is required to satisfy 


Mfj - ^M 2 6ij = -^Sij ( 10 ) 

where Mfj = , M 2 = and S l] is the rate of strain tensor, defined as 

[diUj+djUi]/2. The matrix M can be constructed as follows: in incompressible flow, 
S is a sym metric matrix with eigenvalues that sum to zero; it can be diagonalized 
by the unitary matrix U of eigenvectors: 


S = U diag[A,, A 2 , A 3 ] • e U 


( 11 ) 
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where Ai > A2 > A3 are the eigenvalues in decreasing order. They satisfy 


Ai + A2 + A3 — 0. 


( 12 ) 


In terms of these eigenvalues 

M = U • diag ^0, y / Aj — A2, \/Ai — A3 j • *U . (13) 

By virtue of (12) and (13), M 2 = 3Ai. 

The second moment equation for (8) is derived by applying the rules 


dWi = 0 

dW,dW } = dtSij (14) 

UjdW, = 0 . 

Using these to compute dujuj/dt for (8) shows that (7) is formally the exact second 
moment equation of the stochastic process (8) (Durbin and Speziale 1993). The only 
issue is whether the Langevin equation (8) is well-defined. 

If the coefficients are all bounded, then (8) will be well-defined if the square roots 
are real valued. This requires c a and c 0 to be non-negative, c, is an empirical 
constant of the model (7) that can be chosen to be positive. c 0 is a function given 
by (9). The condition that co > 0 translates to 

P 3 

ci > 1 — (02 + C3) h — c,Ai (15) 

£ £ 

where A/ 2 = 3Ai has been used and the time-scale has been set to T = k/e. The 
right side of this inequality depends on the flow. As long as (15) is met, solutions to 
the general linear model are certain to be realizable; however, (15) is not guaranteed 
to be met, so the general linear model need not always have realizable solutions. 
Our experience with models that are special cases of (7) is that (15) is usually 
satisfied. 

In extreme cases, far from equilibrium, unrealizable solutions to the general linear 
model can be generated. However, the present analysis shows that one easily can 
modify the model to guarantee realizable solutions: just alter c\ so that (15) is 
always satisfied. For instance, if C\ = 1.8 is usually a satisfactory value, then 
replacing this constant by the function 


c\ 


= max 


P 3 

1.8, 1 — (c2 + C3)— + 


(16) 


will enforce the inequality (15). This ensures realizability in extreme cases without 
altering the model in usual cases. 
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FIGURE 3. Solution trajectories of the LRR model projected onto the second- 

in wi ant- third- invariant plane. • • •, original LRR model; , modified, realizable 

LRR model. 


Launder et ai (1975) imposed certain symmetry and normalization constraints 
that led them to specialize (7) by setting 


C-2 


c + 8 _ 8c - 2 

~T r ; C3 = ~TT 


60c — 4 
55 


where c is a constant they set to 0.4 — this is the LRR model, 
condition (16) becomes 


(17) 

With (17), the 


c i = max 


1.8, 1 -0.873- +0.545A! 


(IS) 


Figure 3 shows trajectories of the LRR model for homogeneously sheared turbu- 
lence, with the initial condition 


Sq 2 /e = 20; b u = -0.27; b 22 = -0.33; b 33 = 0.6; b 12 = 6 13 = b 23 = 0. 


These values give P/e = 0, -II = 0.27 and III = 0.053, where II and III are the 
invariants of the Reynolds stress tensor. The trajectories in figure 3 are solutions to 
the LRR model projected into the second- third invariant plane. They start in the 
upper right corner of the triangle and ultimately are attracted to the equilibrium 
point in the lower center or the triangle. The initial value of 1 - 0.873P/; + 0.545Ai 
is 3.73 (note that Aj = Sq 2 /As) so that (15) is violated if cj = 1.8. Realizable 
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solutions must remain inside the curvilinear triangle of figure 3. The dotted curve 
shows that the LRR model exits the realizable region for this initial condition. Thus, 
although (15) is a sufficient but not necessary condition, it provides an insight into 
the existence of unrealizable trajectories. The solution trajectory remains inside 
the realizable region when the modification (18) is applied, as shown by the dashed 
curve in figure 3. 

3. Future plans 

Development of numerical tools for solving Reynolds-stress models in complex 
geometries is under way. These are needed so that we can test and further develop 
models. The computer program being developed is based on the incompressible 
Navier-Stokes solver (INS-2D) written by Stuart Rogers at NASA Ames. This is a 
general geometry code and the Reynolds stress solver is being written so that models 
can be implemented readily. The program under development is an extension of 
that used to calculate figure 1. 

The greatest potential utility of Reynolds stress transport models is in strongly 
non-equilibrium flows. However, certain difficulties remain. For instance, all exist- 
ing models predict incorrect rates of relaxation toward equilibrium in highly per- 
turbed flows. Future work will include this area. 

The analysis of realizability raises some intriguing questions. The realizability of 
models for non- homogeneous flows is largely unexplored. Although I have consid- 
ered this issue for the type of near-wall models that I have developed, I have no 
mathematical understanding of why they seem to produce realizable solutions. The 
ellipticity of non-local wall effects suggests using maximum principles. 

Stochastic differential equations are popular models for Lagrangian dispersion 
calculations. The present use of stochastic equations was solely to analyze statistical 
moment models. This work can be made the basis of self-consistent Reynolds-stress 
and Lagrangian dispersion models. 
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